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This paper begins the process of verifying and validating computational fluid dynam- 
ics (CFD) codes for supersonic retropropulsive flows. Four CFD codes (DPLR, FUN3D, 
OVERFLOW, and US3D) are used to perform various numerical and physical modeling 
studies toward the goal of comparing predictions with a wind tunnel experiment specifi- 
cally designed to support CFD validation. Numerical studies run the gamut in rigor from 
code-to-code comparisons to observed order-of-accuracy tests. Results indicate that this 
complex flowfield, involving time-dependent shocks and vortex shedding, design order of 
accuracy is not clearly evident. Also explored is the extent of physical modeling necessary 
to predict the salient flowfield features found in high-speed Schlieren images and surface 
pressure measurements taken during the validation experiment. Physical modeling studies 
include geometric items such as wind tunnel wall and sting mount interference, as well 
as turbulence modeling that ranges from a RANS (Reynolds- Averaged Navier-Stokes) 2- 
equation model to DES (Detached Eddy Simulation) models. These studies indicate that 
tunnel wall interference is minimal for the cases investigated; model mounting hardware 
effects are confined to the aft end of the model; and sparse grid resolution and turbulence 
modeling can damp or entirely dissipate the unsteadiness of this self-excited flow. 


Nomenclature 


a 

Angle of attack 

m 

Slope 

At 

Time step 

N 

Number of grid nodes 

m 

mass flow rate 

PO.jet. 

Jet stagnation pressure 

e 

Angle measured from top of model 

Poo 

Freestream static pressure 


anti-clockwise when viewed along the 

R 

Model radius, 2.5" 


positive x direction 

r 

Radius from x axis 

C A 

Coefficient of axial force 

R 2 

Coefficient of determination 

Cp 

Coefficient of pressure 

Re 

Reynolds number 

c D 

Coefficient of drag 

Rljet 

Jet stagnation temperature 

C L 

Coefficient of lift 

Too 

Freestream static temperature 

C t 

Coefficient of thrust 

X 

Distance from spherical nose cap 

D 

Model diameter, 5" 


along model axis 

f 

Frequency 

y 

Distance from model axis out the 

h 

Mesh spacing 


starboard side 

L 

Model extent in rr-direction, 10.5" 

z 

Distance from model axis upward 
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I. Introduction 


Future Mars missions will demand entry, descent, and landing (EDL) performance beyond the capabili- 
ties of current parachute-based deceleration systems due to much higher payload masses required for human 
exploration. One alternate deceleration technique is the use of propulsive deceleration during the supersonic 
phase of the entry trajectory as opposed to merely the terminal-landing phase. This supersonic retropropul- 
sion (SRP), or the use of engine thrust directed into the oncoming supersonic freestream flow, has become 
a focus of the EDL community. 1-5 

Because Earth-based testing facilities cannot simultaneously simulate all Mars flight parameters, a portion 
of the Mars entry aerodynamic database will have to come from CFD prediction. The nuclear power and 
weapons communities have long faced a similar circumstance: a lack of, or outright ban on, experimental 
testing. They have spent decades researching the verification, validation, and uncertainty quantification 
steps necessary to build confidence in their numerical predictions. 6 

This trend towards flight vehicle design based solely on CFD is discussed by Aesclrliman and Oberkampf, 7 
who also describe a systematic approach for designing CFD validation experiments. To date, only a handful 
of ground-based retropropulsion experiments have been conducted, 8-10 and they were not designed for CFD 
validation as revealed by recent attempts to simulate them with CFD. 11-15 Aeschliman and Oberkampf 
warned of these difficulties: 

We consider unsatisfactory the common practice of attempting to validate codes using published 
data obtained for some purpose unrelated to code validation. Almost inevitably, critical infor- 
mation required by the code, boundary and initial conditions especially, will be unavailable. 

To rectify this, a CFD validation experiment has been designed, 16 phase I testing has been conducted in the 
NASA Langley 4'x4' Unitary Plan Wind Tunnel (UPWT), 17 and phase II, an entry in the NASA Ames 9'x7' 
Supersonic Wind Tunnel, is scheduled to be complete by the end of 2011. While the rigorous uncertainty 
quantification analysis has yet to be completed for phase I, initial results are available for comparison and 
will be revealed toward the end of this paper. The bulk of this paper, however, chronicles many of the 
numerical verification and physical modeling studies done to prepare for predicting the results from this 
CFD validation experiment and goes beyond those documented during the design of the CFD validation 
experiment itself. 16,18 

The following section contains a brief description of the validation experiment and the case chosen as 
the initial focus for comparison. Next, each of the four CFD codes employed, Dplr, 19 Fun3D, 20,21 Over- 
flow, 22 and Us3D, 23 are discussed along with code-specific verification, solution verification, and numerical 
sensitivity studies. Afterward, several modeling effects are considered: wind tunnel wall interference, sting 
effects, and turbulence modeling. Finally, CFD predictions are assessed against the preliminary experimental 
data and followed by concluding remarks. 
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II. Validation Experiment Case Considered 


The 5"-diameter, 70-degree sphere-cone-cylinder wind tunnel model as mounted in test section 2 of the 
NASA Langley UPWT, 24 is shown in Figure 1 while the surface pressure instrumentation layout is shown in 
Figure 2. Thanks to an interchangeable nozzle-plug design, the model can have a variety of configurations. 



Figure 1: Mounting apparatus and the 5"-diameter model with 3 nozzles as installed in test section 2 of the NASA Langley 
Unitary Plan Wind Tunnel. 24 


During the phase I test campaign, 0-, 1-, 3-, and 4-jet configurations were tested at Mach 2.4, 3.5, and 4.6 
with angles of attack ranging from —8° to 20° and Reynolds number per foot of 1.5M for Mach 4.6 and 1M 
for Machs 2.4 and 3.5. Nominal jet thrust coefficients, Ct, spanned 0.5 to 3 with isolated runs as high as 6. 
For more details of the experiment, please refer to the companion paper by Berry et al. 17 

The experimental run chosen for comparison in this paper, Run 165, represents the simplest jet-on case 
available in that it has a quasi-periodic flowfield structure in which a ring vortex at the tripoint of the 
jet-plume barrel shock, a shear layer formed by large entrained recirculation region, and the jet termination 
shock is shed periodically. The zero angle-of-attack point of this run has the additional, simplifying advantage 
of being largely axisymmetric and almost perfectly periodic with a frequency around 2 kHz. 

Run 165 freestream conditions are Mach 4.6, 1.5M Re/ft, and 65 K air with a nozzle pressure ratio, 
Po,jet / Poo , of 7,724 and a nozzle temperature ratio, Tojet/Toc, of 5.34. For the 4:1 area ratio conical nozzles, 
these conditions produce a thrust coefficient, Ct, of 1.97 and a mass flow rate, to, of 0.62 lbm/s. 

Figure 3 on page 5 shows tableaux of high-speed Schlieren movie still frames for 0° and 12° angle of 
attack. If viewed in electronic form with a viewer that supports embedded Adobe Shockwave™ Flash 
movies, selecting a tableau will play the corresponding movie from which the stills were extracted. To exit 
a movie, select it and press the escape key. 
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180 ° 


(a) Forebody as viewed from the front of the model. 
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(b) Aftbody shell as viewed from the port (— y) side of the model. 

Figure 2: Model surface instrumentation layout. The small, numbered open circles (o) are pressure taps while the small, 
numbered filled circles (•) are 40 kHz pressure transducers. The circular outlines in (a) are the locations of the nozzle 
inserts, which provide options for zero to four jets. 


126 


0 129 

0 130 
^13T 

o' 


147 


132 


-.1 33 


144 

o'°°o 

134 


135 


-©- 


157 


-St- 


148 


-Si- 


158 


-©- 


149 


-©- 


154 

m — 


159 


- 0 - 


164 

— Qj 


150 


-©- 


151 

-©— 


160 


-©- 


-Gh— 


168 

-G 


169 

- 0 — 


170 

-79 — 


174 


707 


136 


o 
o' 

q 138 o 5 


137 


155 


-e- 


161 


- 0 - 


162 


165 


-©- 


171 

-79 


172 

-79 — 


175 


-0— 


.127 


139 

140 


152 

- 0 — 


153 

-0 


141 


o 

0142 
143 146 

coo 


156 


-0- 


163 


-St- 


166 


-St 1 


173 

- 70 — 


176 


- 0 - 


4 of 31 


American Institute of Aeronautics and Astronautics 




(b) a = 12°. 

Figure 3: Schlieren movies for Run 165 from Ref. 17. The static image shows a tableau of four representative still frames. 
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III. Numerical Methods, Sensitivities, and Verification 

Four types of Reynolds-Averaged Navier-Stokes (RANS) CFD codes were used to simulate SRP flows: 
cell-centered overset structured grid (Dplr), node-centered unstructured grid (Fun3D), node-centered over- 
set structured grid (Overflow), and cell-centered unstructured grid (Us3D). All codes simulated time- 
dependent, fully-turbulent flow. The following subsections contain a brief summary of each code, how each 
was run, and a description of code and solution verification 6, 25 studies performed to date. 

A. DPLR 

The Dplr (Data Parallel Line Relaxation) CFD code 19 
is a parallel, structured multi-block, finite-volume code 
with overset grid capability that solves the RANS equa- 
tions for continuum flow, including finite-rate chemistry 
and thermal non-equilibrium. In the present study, the 
thermally- and calorically-perfect RANS equations for 
air are solved implicitly with lst-order time accuracy. In- 
viscid fluxes are formed via a modified Steger- Warming 
flux vector splitting 26 with 3rd-order Monotone Upwind 
Scheme for Conservation Laws (MUSCL) extrapolation 27 
subject to a minmod limiter 28 and 2nd-order trapezoidal 
rule flux integration. The viscous fluxes are computed 
with 2nd-order spatial accuracy with a central difference 
approach. For the present analysis, the Menter’s Shear- 
Stress Transport (SST) turbulence model was employed 
with a vorticity-based production term 29 and no com- 
pressibility corrections. 

A design-order verification study was conducted for 
an inviscid, axisymmetric nozzle problem described by 
Roy. 30 This problem is smooth in the sense that it does 
not contain shocks, which reduce the scheme’s design or- 
der of accuracy per Godunov’s theorem. Figure 4 shows 
a plot of L2 error versus grid refinement factor, h , for 
2nd-order and 3rd-order MUSCL flux reconstructions. 

Both options demonstrate nominally 2nd-order conver- 
gence for this smooth flow with 1.9 and 1.7 slopes, re- 
spectively. The 3rd-order scheme has slightly lower over- 
all error but begins to degrade as the mesh is refined. 

3rd-order reconstruction was used for subsequent results 
shown in this paper. 

A topology and grid sensitivity study was performed for Run 165, a = 0° conditions with three grids: 
Grid 1 with 84M nodes, Grid 2 with 18M nodes, and Grid 3 with 53M nodes. Grids 1 and 2 differ with respect 
to the topology of the overset regions in the jet interaction region while Grid 3 is a uniform refinement of 
Grid 2. Grid 2 provided better resolution of the plume region than Grid 1 with only 20% of the points. Based 
on comparison of the two solutions, this adjustment demonstrated that points far from areas of interest were 
not influential. 

Figure 5 shows a comparison of predicted surface pressures by all three grids along the model centerline 
(0 = 0°, 180°). As the overset topology is changed and the grid is refined, the solutions show improvement 
with respect to the amount of asymmetry predicted for this symmetric flow as indicated by the decreasing 
spread between the upper and lower centerline values. Grid 3 was considered the best compromise and was 
used for the predictions in the Preliminary Comparisons section on page 20. 



Figure 4: Dplr observed order-of-accuracy study for 
axisymmetric nozzle flow. 
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0.04 



r/R 

(a) Forebody only. 



x/D 

(b) Entire body. 

Figure 5: Dplr. grid resolution and topology study. 
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Results from a time step resolution study are shown in Figures 6 and 7 where Figure 6 is a detailed view 
of first portion of Figure 7. Both figures show simulated pressure traces of two forebody tap locations, a 105 
and 106, for a period of time with up to four different time steps: the nominal 0.32 /is, one ten times larger, 
another ten times smaller, and in the case of Figure 6, another a hundred times smaller. As shown by both 
figures, the time step ten times larger than the nominal time step shows immediate and significant departure 
from all others while the difference between two smallest time steps, a tenth and a hundredth as small as 
the nominal (available only in Figure 6) is imperceptible on these expanded C p -scale plots. Figure 7’s five 
times longer time duration and much larger C p range that encompasses typical fluctuations shows that the 
difference between the nominal time step and one ten times smaller is negligible, and so the nominal time 
step was used for the rest of the study. 





Figure 6: Dplr time step resolution study: short duration Figure 7: Dplr time step resolution study: longer duration 
with time steps spanning four orders of magnitude. with time steps spanning three orders of magnitude. 


’In the full study, ten pressure tap locations were examined and all showed behavior similar to those shown here. 
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B. FUN3D 


Table 1: Fun 3D spatial order study: 
number of mesh points and observed er- 


The Fun3D (Fully Unstructured Navier-Stokes Three-Dimensional) suite of codes contains a node-based 
finite- volume flow solver. 20 For this study, Fun3D was run with a selectively-dissipative version of the Low- 
Dissipation Flux Splitting Scheme (LDFSS) inviscid flux function 31 and a modified Van Albada limiter 32 
according to Vatsa and White. 33 The compressible RANS equations are loosely coupled to Menter’s SST 
turbulence model with a vorticity-based production term 29 or a modified Delayed Detached Eddy Simulation 
model that uses the Spalart-Allmaras (SA) near-wall model, 34 which will subsequently be referred to as SA- 
DES. Node-based conservative variables are computed by driving a 2nd-order accurate spatial residual to 
steady-state with a point-implicit iterative method. A modified, optimum 2nd-order backward difference 
formula (BDF) scheme is used in conjunction with a temporal error controller that assures design order. 35 

In addition to previous code verification studies and their as- 
sociated automatic regression test suites, 36 solution verification in 
the form of observed spatial order of accuracy was attempted by 
monitoring two error quantities as a function of grid resolution for 
unsteady Run 165, a = 0° conditions. All grids are of the same 
family in that the spacing field used by the grid generator was uni- 
formly refined from coarse to fine grids. The viscous spacing and 
stretching rates normal to the wall, however, were held constant. 

This allows more surface grid resolution as the grids are refined, 
but the wall-normal distribution remains fixed. Table 1 provides 
the range of grid sizes considered for this study. Also tabulated 
are the two error measures used to investigate the spatial order 
behavior. The first, coefficient of lift (Cl) should be zero for this 
zero angle of attack case and the second is a derived error by com- 
puting the difference between the coefficient of drag (Cjj) on the 

finest and the current grid. For this unsteady problem, both of these measures are inherently noisy due to 
the time-averaging necessary to determine these quantities. 

These two errors are plotted as a function of the number of nodes to the two-thirds power, N~ 2 ^, in 
Figure 8. On these plots, a 2nd-order scheme (the smooth-flow design order for Fun 3D) has a slope of 
one. b As anticipated, both subfigures show very noisy results for this non-smooth flow with such a loose 
family of grids. The lift error behavior is especially bad for finer grids. The non-monotonic slopes may be 
due to different feature sets becoming resolved as the grid is refined as shown by the simulated Schlierens in 
Figure 9 or simply that this complex, unsteady flow is beyond such crude measurement techniques. Perhaps 
running even finer grid levels would more clearly indicate a transition to lst-order, which is hinted at in 
subfigure 8b and is expected for non-smooth flow as the grid is refined. 30 

During this investigation, a simulated Schlieren capability 0 was added to Fun 3D and provided a revealing 
visual depiction of grid resolution effects — see Figure 9, which shows simulated Schlieren images for three of 
the grids used in the spatial order study. As the grid is refined, shocks and shear layers sharpen. Due to 
lack of resources, however, only the medium-resolution grid (28M nodes) was used for the final simulations. 

For Fun 3D runs, the solution was advanced in time with an optimal 2nd-order modified BDF scheme 35 
with a non-dimensional time step of 0.003 at a CFL (Courant-Friedrichs-Lewy) number of five. This time 
step, which gave approximately 200 steps per cycle, corresponds to a physical time step of 2.36 /zs. A 
maximum of 20 sub-iterations were used for each time step. A temporal error controller was used to assure 
temporal accuracy. If the ^-momentum and turbulent kinetic energy L2 error norms dropped an order of 
magnitude below the temporal error estimate, the solution was advanced to the next time step. 


Nodes 

C L 

C D C-D, finest 

3.8M 

2.6 x 10“ 3 

8.78 x 10' 4 

6.7M 

1.4 x 10~ 3 

6.52 xlO' 4 

10. 8M 

7.1 x 10~ 4 

8.72 xlO' 5 

19. 7M 

1.6 x 10~ 4 

1.17 x 10~ 4 

28. 4M 

3.6 xlO' 4 

8.65 xlO' 5 

43.4M 

4.8 x 10' 4 

8.93x10'® 

72. 6M 

1.8 x 10“ 5 

- 


b The N 2 / 3 combines a measure of the mesh spacing, h ~ N 1 / 3 , and 2nd-order behavior with respect to spacing, h 2 . 
c Simulated Schlieren images are generated via volume integration as described by Yates. 37 
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-3 


-2 



( a ) Cl- (b) Cd Cl), finest- 

Figure 8: Fun3D observed spatial order behavior. On this plot, unity slope represents 2nd-order 
convergence. 



(a) 3.8M nodes. 


(b) 28.4M nodes. 


(c) 72. 6M nodes. 


Figure 9: Visualizing grid resolution effects via 700 x 700-pixel simulated Schlieren (Fun3D). 
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C. OVERFLOW 


Overflow 2 (OVERset grid FLOW solver) is an implicit RANS 
flow solver that utilizes structured overset grids. 22 For the cur- 
rent work, the HLLE++ numerical flux function 38 with the Van 
Albada limiter 32 was used for spatial terms, and the Symmetric 
Successive Over Relaxation (SSOR) algorithm with Newton sub- 
iterations for temporal terms. All viscous terms were included, 
and turbulence was modeled with either Menter’s original SST 
model 39 subject to Wilcox’s realizability constraint 40 or a SST- 
DES hybrid model. 41 The overall scheme is 2nd-order accurate 
in space and time. The calculation of inviscid fluxes for both the 
flow solver and the turbulence model use 3rd-order accurate MUSCL reconstruction. 

A grid sensitivity study was conducted with the three axisymmetric grids described in Table 2. The grid 
spacing levels are based on a global refinement factor where as the refinement factor decreases, refinement 
increases. A couple views of the finest grid are shown in Figure 10. Run 165, a = 0° conditions were used 
and the domain was modeled without tunnel walls. 


Table 2: Overflow axisymmetric grid 
sensitivity study. 


Name 

Nodes 

Refinement/dir 

coarse 

169,710 

1.00 

medium 

641,127 

1.94 

fine 

2,656,902 

3.96 



(a) Midrange view. (b) Close up. 

Figure 10: Fine grid used for Overflow axisymmetric grid convergence study. 


Figure 11 shows time-averaged solution comparisons for three grid resolutions. For C p on the model face, 
the fine and medium grids are similar while the coarse grid is low. Downstream of the model shoulder shows 
a lot of difference in C p between refinement levels. Meanwhile, the shock standoff distances of all refinements 
are similar with the coarse grid being slightly farther forward. 

To gauge the differences in unsteadiness captured by the three grid resolutions, the frequency of the drag 
force oscillations were computed as 1.7, 2.0, and 2.3 kHz from coarse to find grids, respectively. The fine and 
medium grids are similar in oscillation period and drag amplitude while the coarse grid has both smaller 
period and amplitude. The higher the refinement, the more detail captured; and similarly, the amount 
of unsteadiness increases with grid refinement. For the fine and medium grids, the average C p on frontal 
face is very similar and the average C p on side body oscillate near each other. Shedding from the tripoint 
region increases average C p on model face; and for time-averaged results, the medium grid yields answers 
comparable to the fine grid at a fraction of the cost. 

Another study was conducted to measure spatial order of accuracy. Here, a family of five axisymmetric 
grids were used as shown in Table 3. As with Fun3D, the error in Cd with respect to the finest mesh was 
used to monitor convergence. Figure 12 shows the error as a function of mesh spacing. A linear fit of the 
points gives a slope of 1.0, indicating 2nd-order spatial accuracy. 
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Figure 11: Overflow time-averaged solution comparison for three grid resolutions. 


Table 3: Overflow axisymmetric grids 
for spatial order study. 


Name 

Nodes 

Refinement/dir 

Grid 1 

169,710 

1.00 

Grid 2 

287,046 

1.30 

Grid 3 

641,127 

1.94 

Grid 4 

1,098,819 

2.54 

Grid 5 

2,656,902 

3.96 



Log 10 h 2 


Figure 12: Overflow observed spatial order of accuracy for axisymmetric simulation of Run 165, a = 0°. 
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Next, the number of subiterations necessary during time advancement was investigated. Starting from 
time step 11,500 on the medium grid with a 1.71 /zs time step, 5, 10, and 20 subiterations were run. As 
expected, the drop in subiteration L2 residuals increases with number of subiterations. No variation in drag 
or frequency was seen, and therefore five subiterations were used henceforth. 

Finally a time step sensitivity study was performed where the time step was varied though three orders 
of magnitude: 17.1, 1.71, and 0.171 /zs. Starting at the same point in time, each was run for 2,760 /zs as 
shown in Figure 13. The largest time step, 17.1 /zs, has lower frequency than the smaller time steps, For 
time steps of 1.71 /zs and 0.171 /zs, drag behavior is very similar. Time-averaged shock standoff distances, 
as shown in Figure 14, are not very sensitive to time step. The two smallest time steps are nearly identical 
in time-averaged C p distributions. Drag oscillation frequencies as a function of decreasing time step size 
were 1.17, 1.83, and 1.99 kHz, respectively. As time step is reduced, the frequency of the tripoint oscillation 
increases. In conclusion, to get average flow field properties, a time step of 1.71 /zs is sufficient. 



Figure 13: Overflow Cq as a function of time step. 
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x/D 

Figure 14: Overflow time-averaged shock standoff distance as a function of time step. 


D. US3D 

The Us3D (Unstructured Three-Dimensional) code 23,42 is a parallel, Navier-Stokes, finite- volume solver 
for unstructured grids. It shares roughly the same features, numerical methods and physical models as 
Dplr including the finite-rate chemistry and thermal non-equilibrium with a major difference being the 
unstructured grid capability of Us3D. For a given problem using identical grids and modeling options, the 
two codes produce the exact same solutions to machine precision. 

Us3D is a relatively new code and is still in development phase. As of this writing, two equation 
turbulence models are yet to be implemented. In the current study, SA-DES turbulence modeling is used. 43,44 
Both inviscid and viscous fluxes are evaluated with 2nd-order spatial accuracy, and backward Euler time 
integration provided lst-order temporal accuracy. 

Unfortunately, Us3D grid and time resolution sensitivity studies have not yet been completed for the 
problem studied herein. The grid employed consisted of 31.4 million hexahedral cells, and the time step was 
fixed at 1.5 /is. 
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IV. Physical Modeling Studies 


In addition to code verification and solution verification discussed above, three aspects of physical model- 
ing were explored: wind tunnel wall interference, the effect of the model support hardware and high-pressure 
feedline, and turbulence modeling. These were explored to determine what level of physical modeling was 
necessary to capture the flow physics. 

A. Tunnel Wall Effects 

Figures 15 and 16 show the effects on the time-averaged surface pressures of modeling inviscid tunnel walls 
or not with Overflow. For these angles of attack, there are no substantial differences in surface C p 
and the differences seen in the forebody distribution of Figure 16a is assumed to be due to insufficient 
time-averaging. This 12° angle of attack case contains large-scale unsteadiness of the entire flowfield, and 
the original practice of averaging over several thousand time steps was later found to be insufficient when 
monitoring root-mean-square fluctuations. 




(b) Entire body. 


Figure 15: Effect of inviscid tunnel walls on forebody C p for Run 165 at a = 0°. 




(a) Forebody only. (b) Entire body. 


Figure 16: Effect of inviscid tunnel walls on surface C p for Run 165 at a = 12°. 
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B. Sting Effects 


On the sting, the high pressure feedline and instrumentation wiring bundles caused protuberances close to 
the model base. Since a subsonic wake region existed in this area, it was possible that influence from the 
protuberances traveled upstream and affected pressure tap data. To quantify any discrepancy in pressure 
tap readings due to the sting, an overset mesh was created from a simplified sting geometry, and cases were 
simulated using Overflow. The goal was to show whether including the sting in the CFD models was a 
requirement for good comparison to surface pressure test data. 

The model was greatly simplified to get a rough idea of lst-order effects as shown in Figure 17. Dimensions 
were estimated from photographs; and while not exact, they should be close enough for the purposes of this 
study. As seen in the figure, the sting was modeled only to a certain point downstream. Cavities under the 
white fiberglass tape holding the feedlines to the sting were considered solid, and the instrumentation wiring 
bundle was not modeled. 



Figure 17: Model mount modeling. 


Four cases were simulated at Mach 2.4, Re/ ft = 1 M, and Ct = 3. They were a single nozzle configuration 
at zero angle of attack and a triple nozzle configuration at angles of attack -8, 0, and 12 degrees. Each case 
was run with and without the sting, and surface pressure comparisons were conducted. A sample flowfield 
with the sting is shown in Figure 18 for the tri-nozzle configuration at a = 0°. 



Figure 18: Sting flow. 
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Comparisons of time-averaged surface pressures for the tri-nozzle configuration at zero angle of attack 
case are shown in Figure 19. This figure shows that the sting raised the base pressure on the model. This 
influence traveled upstream far enough to affect the three pressure ports that were the furthest aft on the 
model, taps 174-176 — see Figure 2 on page 4. The compression at the model base due to the sting also 
caused a decrease in drag. This study also showed changes in the unsteady fluctuations that affected the 
instantaneous surface pressures, but the time-averaged values were unaffected. 


Free 



Figure 19: Sting effects. 

Because the sting effects were limited to the base and only three aft pressure ports, the effort of modeling 
the sting had limited returns, and was deemed unnecessary for pressure comparisons over the remainder of 
the body. In future tests, however, where force and moment quantities are measured, the discrepancy in the 
model base pressure will be an important consideration because the sting causes increased base pressure and 
thus directly affects the model’s drag. 
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C. Turbulence Modeling 


The turbulence models considered span: Spalart- 

Allmaras (SA) 43 and two versions of Menter’s k-ui 
Shear-Stress Transport (SST) 29,39 Reynolds- Averaged 
Navier-Stokes (RANS) models; Spalart’s Delayed 
Detached Eddy Simulation (DDES) model; 45 and 
Sanchez-Rocha et al’s Hybrid RANS-Large Eddy Sim- 
ulation (HRLES) model. 46 

Figure 20 shows turbulent eddy viscosity levels pre- 
dicted by Overflow relative to laminar (bulk) vis- 
cosity for three configurations of Menter’s original SST 
turbulence model: (1) the baseline model, 39 (2) the 
baseline model with a compressibility correction, 47 and 
(3) the baseline model with Wilcox’s realizability con- 
straint. 40 The turbulent eddy viscosity levels for the 
baseline model are high between bow and terminal 
shocks, which is a known, non-physical artifact. The 
baseline model with compressibility correction is simi- 
lar in behavior, but has lower levels of eddy viscosity 
and the flow is more unsteady. Eddy viscosity levels for 
the baseline model with a realizability constraint are 
lower than the other models around the shocks, and 
higher in the wake. The realizability constraint elim- 
inates the spurious post-shock eddy viscosity produc- 
tion and better matches wind tunnel plume interaction 
behavior. 

Figure 21 shows a comparison of Fun3D time- 
averaged surface pressures between the SST turbulence 
model and SA-DES along the upper and lower cen- 
terlines for Run 165 at zero angle of attack. In the 
forebody region, and especially near the jet exit at 
r/R = 0.1, SA-DES recovers higher pressures than 
RANS. Based on solution behavior during grid con- 
vergence studies, this increase in pressure recovery is 
thought to be related to the extent that the jet interac- 
tion is dissipated. For the SA-DES solution, the RANS 
submodel is only engaged in the near-wall regions and 
thus does not add turbulent eddy viscosity to the jet 
interaction region, leading to less overall dissipation as 
compared to the pure RANS solution. 



(a) Standard model. 



(b) Compressibility correction. 



(c) Realizability constraint. 


Figure 20: Relative eddy viscosity levels for various 

Menter SST turbulence model variations. 
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(b) Entire body. 

Figure 21: Average surface pressure comparisons between Menter's SST turbulence model and SA-DES for Run 165, 
a = 0° (Fun3D). 
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V. Preliminary Comparisons 


Given that the uncertainty analysis for the validation experiment was not complete at the time of writing, 
this section will concentrate on code-to-code comparisons but will also include some preliminary experimental 
data without uncertainty estimates or corrections for wind tunnel flow angle, bias errors, or random errors. 
Three forms of comparisons will be made in this section: grid spacing distributions used by each code, 
prediction of the salient flowfield features and the experimentally observed frequency of oscillation, and 
finally a comparison of time-averaged surface pressures with raw experimental surface pressure data. 

The first comparison is the grid spacing distributions used by each code. Figure 22 depicts a grid slice 
in the pitch plane focused on the jet interaction region. To be able to compare relative spacings between 
the grids even though Fun3D’s tetrahedral grid is sliced into triangles and quadrilateral faces, the slice is 
colored by an isotropic spacing, h , which is the cube root of the cell volume. Each grid has clustering around 
the jet interaction area with Dplr having the finest with 0.005 diameter spacing, followed by Overflow, 
Us3D, and Fun3D. 



(c) Overflow. (d) Us3D. 


Figure 22: Grid spacing distributions in the symmetry plane where h is defined as (cell volume) 1 ^ 3 . 
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Figure 23 shows the forebody surface grid used by each code, where h is now defined as the square root 
of the cell area. Here the distributions are more varied with the Dplr grid having more resolution around 
the jet exit than other grids but considerably less in the shoulder region. The Fun3D grid appears to have 
considerably more shoulder refinement than the others, but this an artifact of h in that Fun3D’s grid is 
isotropic and the rest have circumferentially-stretched anisotropic cells. On the forebody, Overflow has 
the smallest spacing, followed by Dplr, Fun3D, and Us3D. 



(a) Dplr. 


(b) Fun3D. 




(c) Overflow. 



(d) Us3D. 


Figure 23: Forebody surface grid for each code colored by h, which is (Cell Area) 1 / 2 . 


Figures 24 to 27 on pages 22-25 contain animations meant to simulate the experimental Schlieren movie 
shown in Figure 3a on page 5. Dplr and Us3D show Mach contours in the pitch plane, Fun3D has a 
simulated Schlieren via integrated density deflections through the computational volume similar to Ref. 37, 
and Overflow presents a simulated focused shadowgraph via density gradients in the pitch plane. While all 
four codes capture the periodic vortex ring shedding from the tripoint region, the amplitude of the oscillation 
is considerably smaller for Dplr, which is the only code running the Menter’s SST turbulence model — the 
others were run with DES, which does not inject turbulent eddy viscosity in this region. 
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Run 165. W 
1 - 0 . 000*- 00 s 



Figure 24: Dplr animation of symmetry plane Mach contours and surface C p for Run 165, a = 0°. The static image 
shows a tableau of four representative still frames. 
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Figure 25: F(JN3D-simulated Schlieren and surface C p animation for Run 165, a = 0°. The static image shows a tableau 
of four representative still frames. 
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Figure 26: Overflow animation of symmetry plane density gradients and surface C p for Run 165, a = 0°. The static 
image shows a tableau of four representative still frames. 
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Figure 27: Us3D animation of symmetry plane Mach contours for Run 165, a = 0°. The static image shows a tableau 
of four representative still frames. 
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High-speed pressure transducers were post-processed 
to obtain root-mean-square (RMS)-averaged fluctua- 
tions and a power spectral density function for compar- 
ison with computations. 1 ' For the zero angle of attack 
point of Run 165, the ring vortex shedding frequency, 
which is the dominant oscillation of the jet plume struc- 
ture, is given in Table 4 along with all four codes plus 
an inviscid run by Cart3D. 48 As discussed in Ref. 17, 
depending on where the pressure transducer was located 
on the forebody, the predominant frequency ranged be- 
tween 1.7 kHz and 2.3 kHz. This range appears to cor- 
relate to a non-zero angle-of-attack due to tunnel flow 
angle. The frequency for Dplr and Us3D were ob- 
tained by performing a Fast Fourier Transform (FFT) 
of pressure at several forebody tap locations. Cart3D 
computed the frequency via FFT of the axial force while Overflow and Fun3D measured the steps per 
cycle of the axial force. All codes predict a frequency within the range of the experiment results; and as 
indicated by the inviscid Cart3D result, the ring vortex shedding is captured inviscidly. 

Figure 28 shows the comparison of time-averaged surface pressures between codes and the raw exper- 
imental data from Ref. 17 for Run 165, a = 0°. Only the top and bottom centerline data are shown 
(0 = 0°, 180°). Subfigure 28a shows C p on the forebody as a function of radius from the jet exit at r/R = 0.1 
to the shoulder at r/R = 1 while sub figure 28b spans the entire length of the model as a function of x/D. 
The raw experimental data shows some asymmetry on the forebody and significant asymmetry along the 
side body. As mentioned previously, this is thought to be due to an approximately one degree tunnel flow 
angle. Regardless, the C p levels for this zero angle-of-attack point are quite low as a result of the jet plume 
blocking the incoming flow and enveloping the entire model in subsonic recirculation regions; and so, some 
spread is anticipated on this expanded scale. 

All codes predict similar trends for the forebody pressure distribution, however, the amount of pressure 
recovered on the forebody appears to correlate to grid resolution. As shown in Figure 22 on page 20 and 
Figure 23 on page 21, Overflow has the greatest concentration of points on and around the forebody and 
also predicts the highest pressure recovery on the forebody. Despite having relatively fine spacing in the 
forebody area compared to the grids used for Fun3D and Us3D, Dplr has a diminished peak near the jet 
exit. This is attributed to the SST turbulence model dissipating the periodically-shed ring vortex as shown 
in Figure 24 on page 22. Both Fun3D and Us3D ran SA-DES and agree well on the cone flank. Fun3D 
has a lower peak near the jet exit while Us3D has a lower peak at the shoulder, which correspond to where 
each has fewer grid points compared to the other. 

The side-body pressure levels are similar for all codes run with DES while Dplr is slightly above the rest. 
Again, this may correspond to use of Menter’s SST turbulence model. Note: The Us3D surface pressure 
does not drop near the base because the computational domain ends at the base plane. 

Figure 29 on page 28 shows a similar comparison for the 12° angle-of-attack point. Due to angle of attack, 
the pressures are markedly higher than those for the zero degree case and the deviations in the experimental 
data along a given upper or lower cut are not as pronounced. Similar to the zero degree angle-of-attack point 
and consistent with Figure 21 on page 19, the DES solutions (Fun3D, Overflow, and Us3D) predict higher 
pressure recovery on the lower flank of the forebody than the pure RANS solution (Dplr). 


Table 4: Dominate flowfield frequency comparison. 


Source 

At, s 

Steps/Cycle 

f, kHz 

Experiment 1 

2.5 x 10~ 5 

20 

1.7-2. 3 

Dplr 

3.2 x 10 -7 

1,500 

2.1 

Fun3D 

2.4x 10 -6 

206 

2.1 

Overflow 

1.7 x 10~ 6 

286 

2.1 

Us3D 

1.5 x 10~ 6 

392 

1.7 

CartSD* 

1.2 x 10~ 5 

38 

2.1 


t From Berry et al . 17 
t From Bakhtian . 48 
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Figure 28: Time-averaged surface pressure comparisons for Run 165, a = 0 ° . 
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(b) Entire body. 

Figure 29: Time-averaged surface pressure comparisons for Run 165, oi = 12°. 
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VI. Concluding Remarks 


In this paper we continued the process of verifying and validating computational fluid dynamics (CFD) 
codes for supersonic retropropulsive flows. In this installment, four CFD codes documented various numerical 
and physical modeling studies and made an initial comparisons with each other and with preliminary results 
from a wind tunnel experiment specifically designed to support CFD validation. Numerical studies varied 
in rigor from code-to-code comparisons to sensitivity studies to observed order-of-accuracy tests. Results 
indicate that for this complex, non-smootlr flow, design order was not clearly observed. In addition, various 
physical modeling studies indicated that tunnel wall interference is minimal for the cases considered; the 
model mounting hardware effects are confined to the aft region of the model; and the RANS turbulence 
models damp or entirely dissipate the unsteadiness of this self-excited flow unless sufficient grid resolution 
and realizability constraints are used. To further resolve code-to-code differences, future work should consider 
running the same grid for multiple codes. 
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